﻿using System;
using System.Windows.Forms;
using System.Security.Cryptography;
using System.Globalization;
using System.IO;

namespace TestModCholeskyDecomposition
{
    //Enumerator for describing selected RBF function
    public enum RBFFunction : int
    { EXP = 0, CP_C2 = 1, TPS = 2, MQ = 3 };

    public class Form1 : Form
    {
        /// <summary>
        /// Required designer variable.
        /// </summary>
        private System.ComponentModel.IContainer components = null;

        /// <summary>
        /// Clean up any resources being used.
        /// </summary>
        /// <param name="disposing">true if managed resources should be disposed; otherwise, false.</param>
        protected override void Dispose(bool disposing)
        {
            if (disposing && (components != null))
            {
                components.Dispose();
            }
            base.Dispose(disposing);
        }

        #region Windows Form Designer generated code

        /// <summary>
        /// Required method for Designer support - do not modify
        /// the contents of this method with the code editor.
        /// </summary>
        private void InitializeComponent()
        {
            this.oneTestButton = new System.Windows.Forms.Button();
            this.testsButton = new System.Windows.Forms.Button();
            this.exitButton = new System.Windows.Forms.Button();
            this.label1 = new System.Windows.Forms.Label();
            this.label2 = new System.Windows.Forms.Label();
            this.label3 = new System.Windows.Forms.Label();
            this.nTextBox = new System.Windows.Forms.TextBox();
            this.mTextBox = new System.Windows.Forms.TextBox();
            this.comboBox1 = new System.Windows.Forms.ComboBox();
            this.SuspendLayout();
            // 
            // oneTestButton
            // 
            this.oneTestButton.Location = new System.Drawing.Point(12, 101);
            this.oneTestButton.Name = "oneTestButton";
            this.oneTestButton.Size = new System.Drawing.Size(75, 38);
            this.oneTestButton.TabIndex = 0;
            this.oneTestButton.Text = "One calculation";
            this.oneTestButton.UseVisualStyleBackColor = true;
            this.oneTestButton.Click += new System.EventHandler(this.oneTestButton_Click);
            // 
            // testsButton
            // 
            this.testsButton.Location = new System.Drawing.Point(127, 101);
            this.testsButton.Name = "testsButton";
            this.testsButton.Size = new System.Drawing.Size(75, 38);
            this.testsButton.TabIndex = 1;
            this.testsButton.Text = "Perform tests";
            this.testsButton.UseVisualStyleBackColor = true;
            this.testsButton.Click += new System.EventHandler(this.testsButton_Click);
            // 
            // exitButton
            // 
            this.exitButton.Location = new System.Drawing.Point(245, 101);
            this.exitButton.Name = "exitButton";
            this.exitButton.Size = new System.Drawing.Size(75, 38);
            this.exitButton.TabIndex = 2;
            this.exitButton.Text = "Exit";
            this.exitButton.UseVisualStyleBackColor = true;
            this.exitButton.Click += new System.EventHandler(this.exitButton_Click);
            // 
            // label1
            // 
            this.label1.AutoSize = true;
            this.label1.Location = new System.Drawing.Point(10, 18);
            this.label1.Name = "label1";
            this.label1.Size = new System.Drawing.Size(159, 13);
            this.label1.TabIndex = 3;
            this.label1.Text = "Number of interpolation points n:";
            // 
            // label2
            // 
            this.label2.AutoSize = true;
            this.label2.Location = new System.Drawing.Point(10, 44);
            this.label2.Name = "label2";
            this.label2.Size = new System.Drawing.Size(204, 13);
            this.label2.TabIndex = 4;
            this.label2.Text = "Dimension of the interpolation points l_{x}:";
            // 
            // label3
            // 
            this.label3.AutoSize = true;
            this.label3.Location = new System.Drawing.Point(10, 70);
            this.label3.Name = "label3";
            this.label3.Size = new System.Drawing.Size(112, 13);
            this.label3.TabIndex = 5;
            this.label3.Text = "Radial Basis Function:";
            // 
            // nTextBox
            // 
            this.nTextBox.Location = new System.Drawing.Point(220, 15);
            this.nTextBox.Name = "nTextBox";
            this.nTextBox.Size = new System.Drawing.Size(100, 20);
            this.nTextBox.TabIndex = 6;
            this.nTextBox.Text = "15";
            // 
            // mTextBox
            // 
            this.mTextBox.Location = new System.Drawing.Point(220, 41);
            this.mTextBox.Name = "mTextBox";
            this.mTextBox.Size = new System.Drawing.Size(100, 20);
            this.mTextBox.TabIndex = 7;
            this.mTextBox.Text = "3";
            // 
            // comboBox1
            // 
            this.comboBox1.FormattingEnabled = true;
            this.comboBox1.Items.AddRange(new object[] {
            "EXP",
            "CP_C2",
            "TPS",
            "MQ"});
            this.comboBox1.Location = new System.Drawing.Point(220, 67);
            this.comboBox1.Name = "comboBox1";
            this.comboBox1.Size = new System.Drawing.Size(100, 21);
            this.comboBox1.TabIndex = 8;
            this.comboBox1.SelectedIndex = 0;
            // 
            // Form1
            // 
            this.AutoScaleDimensions = new System.Drawing.SizeF(6F, 13F);
            this.AutoScaleMode = System.Windows.Forms.AutoScaleMode.Font;
            this.ClientSize = new System.Drawing.Size(326, 143);
            this.ControlBox = false;
            this.Controls.Add(this.comboBox1);
            this.Controls.Add(this.mTextBox);
            this.Controls.Add(this.nTextBox);
            this.Controls.Add(this.label3);
            this.Controls.Add(this.label2);
            this.Controls.Add(this.label1);
            this.Controls.Add(this.exitButton);
            this.Controls.Add(this.testsButton);
            this.Controls.Add(this.oneTestButton);
            this.FormBorderStyle = System.Windows.Forms.FormBorderStyle.FixedSingle;
            this.MaximizeBox = false;
            this.Name = "Form1";
            this.ShowIcon = false;
            this.Text = "Test Tool for the Modified Cholesky Decomposition";
            this.ResumeLayout(false);
            this.PerformLayout();

        }

        #endregion

        #region Windows Form objects
        private System.Windows.Forms.Button oneTestButton;
        private System.Windows.Forms.Button testsButton;
        private System.Windows.Forms.Button exitButton;
        private System.Windows.Forms.Label label1;
        private System.Windows.Forms.Label label2;
        private System.Windows.Forms.Label label3;
        private System.Windows.Forms.TextBox nTextBox;
        private System.Windows.Forms.TextBox mTextBox;
        private System.Windows.Forms.ComboBox comboBox1;
        private bool validInput = false;
        #endregion

        //Constructor
        public Form1()
        {
            InitializeComponent();
            nTextBox.KeyDown += new KeyEventHandler(integer_KeyDown);
            mTextBox.KeyDown += new KeyEventHandler(integer_KeyDown);
            nTextBox.KeyPress += new KeyPressEventHandler(integer_KeyPress); ;
            mTextBox.KeyPress += new KeyPressEventHandler(integer_KeyPress);
        }

        #region Forms's methods
        //methods for allowing only integer input
        private void integer_KeyDown(object sender, KeyEventArgs e)
        {
            validInput = true;
            if (
            (e.KeyCode >= Keys.D0 && e.KeyCode <= Keys.D9) ||
            (e.KeyCode >= Keys.NumPad0 && e.KeyCode <= Keys.NumPad9) || e.KeyCode == Keys.Back)
            { validInput = false; }
        }
        private void integer_KeyPress(object sender, KeyPressEventArgs e)
        {
            if (validInput == true)
            { e.Handled = true; }
        }

        //Exit Button
        private void exitButton_Click(object sender, EventArgs e)
        { this.Close(); }

        //Button for performing only the defined test with the input fields
        private void oneTestButton_Click(object sender, EventArgs e)
        {
            int n = 15, m = 3, i;
            RBFFunction function = RBFFunction.EXP;
            try
            { n = int.Parse(nTextBox.Text); }
            catch { n = 15; }
            try { m = int.Parse(mTextBox.Text); }
            catch { m = 3; }
            try
            { function = (RBFFunction)comboBox1.SelectedIndex; }
            catch { function = RBFFunction.EXP; }
            double[][] points = CreatePoints(n, m);
            double[] values = new double[n];
            for (i = 0; i < n; i++)
            { values[i] = Function(points[i], m); }
            string message = PerformLUTest(points, values, n, m, function, true);
            message += PerformModifiedCholeskyTest(points, values, n, m, function, true);
            MessageBox.Show(message);
        }

        //Button for performing all test described in the paper
        private void testsButton_Click(object sender, EventArgs e)
        {
            int i, j;

            //tests parameters
            int[] ns = new int[] { 15, 500, 1000, 1500 };
            int[] ms = new int[] { 3, 5, 10 };
            RBFFunction[] functions = new RBFFunction[] { RBFFunction.EXP, RBFFunction.CP_C2, RBFFunction.TPS, RBFFunction.MQ };
            int nTest = 10;         //Number of repetition of test for each test parameters

            //file for the report
            string file;
            SaveFileDialog saveDialog = new SaveFileDialog();
            saveDialog.Filter = "Text File (*.txt)|*.txt";
            saveDialog.FilterIndex = 1;
            if (saveDialog.ShowDialog() == DialogResult.OK)
            {
                file = saveDialog.FileName;
                if (!String.IsNullOrEmpty(file))
                {
                    TextWriter reporter = new StreamWriter(file, false);
                    reporter.WriteLine("Times are in processor ticks. To calculate the time in milliseconds, you should divide by: "
                        + TimeSpan.TicksPerMillisecond.ToString(CultureInfo.InvariantCulture));
                    reporter.WriteLine("Method\tn\tm\tRBFfunction\tFilling time\tSolving time\tError");
                    //Arrays for the random points for tests
                    double[][] points;
                    double[] values;
                    foreach (int n in ns)
                    {
                        foreach (int m in ms)
                        {
                            for (i = 0; i < nTest; i++)// 
                            {
                                points = CreatePoints(n, m);
                                values = new double[n];
                                for (j = 0; j < n; j++)
                                { values[j] = Function(points[j], m); }
                                foreach (RBFFunction function in functions)
                                {
                                    reporter.WriteLine(PerformLUTest(points, values, n, m, function, false));
                                    reporter.WriteLine(PerformModifiedCholeskyTest(points, values, n, m, function, false));
                                    reporter.Flush();
                                }
                            }
                        }
                    }
                }
                else
                { MessageBox.Show("You did not specify the file. The test will not perform!"); }
            }
        }
        #endregion

        #region Methods for performing the tests
        /// <summary>
        /// Method to perform single test of the LU decomposition.
        /// </summary>
        /// <param name="points">Array of the interpolation points.</param>
        /// <param name="values">Array of the values in the interpolation points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <param name="function">Chosen RBF function.</param>
        /// <param name="messageBox">Boolean to define, if the message is for the messageBox or for the report.</param>
        /// <returns>String with the result of the test.</returns>
        private string PerformLUTest(double[][] points, double[] values, int n, int m, RBFFunction function, bool messageBox)
        {
            #region LU decomposition and solution of the LSE
            double[][] G;
            double[] rhs;
            double[] sol;
            DateTime start = DateTime.Now;
            //filling of the matrix
            switch (function)
            {
                case RBFFunction.EXP:
                    G = FillG_EXP(points, n, m);
                    break;
                case RBFFunction.CP_C2:
                    G = FillG_CP_C2(points, n, m);
                    break;
                case RBFFunction.TPS:
                    G = FillG_TPS(points, n, m);
                    break;
                case RBFFunction.MQ:
                    G = FillG_MQ(points, n, m);
                    break;
                default:
                    G = FillG_EXP(points, n, m);
                    break;
            }
            rhs = FillRHSVector(values, n, m);
            DateTime filled = DateTime.Now;
            sol = LUdecomposition(G, rhs, n + m + 1);
            DateTime solved = DateTime.Now;
            #endregion

            #region Calculation of the error and creating the result string
            //Variables for the error calcualtion
            double lengthError = 0.0;
            double[] errorVector = new double[n + m + 1];
            double[] rhsSol;  //Right hand side of the LSE, calculated with the solution vector and an original matrix G
            //reevaluation of the LSE in case if there were any changes during the decomposition
            rhs = FillRHSVector(values, n, m);
            switch (function)
            {
                case RBFFunction.EXP:
                    G = FillG_EXP(points, n, m);
                    break;
                case RBFFunction.CP_C2:
                    G = FillG_CP_C2(points, n, m);
                    break;
                case RBFFunction.TPS:
                    G = FillG_TPS(points, n, m);
                    break;
                case RBFFunction.MQ:
                    G = FillG_MQ(points, n, m);
                    break;
                default:
                    G = FillG_EXP(points, n, m);
                    break;
            }
            //Calculation of the error of decomposition
            rhsSol = CalculateRHSVector(G, sol, n + m + 1);
            errorVector = CalculateErrorVector(rhsSol, rhs, n + m + 1);
            lengthError = LengthVector(errorVector, n + m + 1);
            string result;
            if (messageBox)     //Report for the messageBox
            {
                result = "LU decomposition test:\nNumber of points n:\t" + n.ToString() + "\tDimension of points m:\t" +
                    m.ToString() + "\tChosen RBF function:\t";
                switch (function)
                {
                    case RBFFunction.EXP:
                        result += "EXP\n";
                        break;
                    case RBFFunction.CP_C2:
                        result += "CP C2\n";
                        break;
                    case RBFFunction.TPS:
                        result += "TPS\n";
                        break;
                    case RBFFunction.MQ:
                        result += "MQ\n";
                        break;
                    default:
                        result += "EXP\n";
                        break;
                }
                result += "Filling ticks:\t" + ((filled - start).Ticks / TimeSpan.TicksPerMillisecond).ToString("f1", CultureInfo.InvariantCulture) +
                    "\tSolving ticks:\t" + ((solved - filled).Ticks / TimeSpan.TicksPerMillisecond).ToString("f1", CultureInfo.InvariantCulture) +
                    "\tError:\t" + lengthError.ToString("f9", CultureInfo.InvariantCulture) + "\n";
                if (n <= 15) //if the number of points is less than 15, the solution and error vectors are added to the report
                {
                    result += "Original RHS Vector: " + VectorToString(rhs, n + m + 1, "f9") + "\n" +
                        "Solved RHS Vector: " + VectorToString(rhsSol, n + m + 1, "f9") + "\n" +
                        "Error RHS Vector: " + VectorToString(errorVector, n + m + 1, "f9") + "\n" +
                        "Solution of the LSE: " + VectorToString(sol, n + m + 1, "f9") + "\n";
                }
            }
            else         //Report for the analysis
            {
                result = "LU:\t" + n.ToString() + "\t" + m.ToString() + "\t";
                switch (function)
                {
                    case RBFFunction.EXP:
                        result += "EXP\t";
                        break;
                    case RBFFunction.CP_C2:
                        result += "CP C2\t";
                        break;
                    case RBFFunction.TPS:
                        result += "TPS\t";
                        break;
                    case RBFFunction.MQ:
                        result += "MQ\t";
                        break;
                    default:
                        result += "EXP\t";
                        break;
                }
                result += ((filled - start).Ticks / TimeSpan.TicksPerMillisecond).ToString("f1", CultureInfo.InvariantCulture) + "\t" +
                    ((solved - filled).Ticks / TimeSpan.TicksPerMillisecond).ToString("f1", CultureInfo.InvariantCulture) +
                    "\t" + lengthError.ToString("f9", CultureInfo.InvariantCulture) + "\n";
            }
            #endregion
            return result;
        }

        /// <summary>
        /// Method to perform single test of the modified Cholesky decomposition.
        /// </summary>
        /// <param name="points">Array of the interpolation points.</param>
        /// <param name="values">Array of the values in the interpolation points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <param name="function">Chosen RBF function.</param>
        /// <param name="messageBox">Boolean to define, if the message is for the messageBox or for the report.</param>
        /// <returns>String with the result of the test.</returns>
        private string PerformModifiedCholeskyTest(double[][] points, double[] values, int n, int m, RBFFunction function, bool messageBox)
        {
            #region Modified Cholesky decomposition and solution of the LSE
            double[][] Phi;
            double[] f;
            double[][] P;
            double[] sol;
            DateTime start = DateTime.Now;
            //filling of the matrix
            switch (function)
            {
                case RBFFunction.EXP:
                    Phi = FillPHI_EXP(points, n, m);
                    break;
                case RBFFunction.CP_C2:
                    Phi = FillPHI_CP_C2(points, n, m);
                    break;
                case RBFFunction.TPS:
                    Phi = FillPHI_TPS(points, n, m);
                    break;
                case RBFFunction.MQ:
                    Phi = FillPHI_MQ(points, n, m);
                    break;
                default:
                    Phi = FillPHI_EXP(points, n, m);
                    break;
            }
            P = FillPMatrix(points, n, m);
            f = Fill_fVector(values, n);
            DateTime filled = DateTime.Now;
            sol = modifiedCholeskyDecomposition(Phi, P, f, n, m);
            DateTime solved = DateTime.Now;
            #endregion

            #region Calculation of the error and creating the result string
            //Variables for the error calcualtion
            double lengthError = 0.0;
            double[] errorVector = new double[n + m + 1];
            double[] rhsSol;  //Right hand side of the LSE, calculated with the solution vector and an original matrix G
            //evaluation of the LSE due to the different structures used in the decomposition
            double[] rhs = FillRHSVector(values, n, m);
            double[][] G;
            switch (function)
            {
                case RBFFunction.EXP:
                    G = FillG_EXP(points, n, m);
                    break;
                case RBFFunction.CP_C2:
                    G = FillG_CP_C2(points, n, m);
                    break;
                case RBFFunction.TPS:
                    G = FillG_TPS(points, n, m);
                    break;
                case RBFFunction.MQ:
                    G = FillG_MQ(points, n, m);
                    break;
                default:
                    G = FillG_EXP(points, n, m);
                    break;
            }
            //Calculation of the error of decomposition
            rhsSol = CalculateRHSVector(G, sol, n + m + 1);
            errorVector = CalculateErrorVector(rhsSol, rhs, n + m + 1);
            lengthError = LengthVector(errorVector, n + m + 1);
            string result;
            if (messageBox)     //Report for the messageBox
            {
                result = "Modified Cholesky decomposition test:\nNumber of points n:\t" + n.ToString() + "\tDimension of points m:\t" +
                    m.ToString() + "\tChosen RBF function:\t";
                switch (function)
                {
                    case RBFFunction.EXP:
                        result += "EXP\n";
                        break;
                    case RBFFunction.CP_C2:
                        result += "CP C2\n";
                        break;
                    case RBFFunction.TPS:
                        result += "TPS\n";
                        break;
                    case RBFFunction.MQ:
                        result += "MQ\n";
                        break;
                    default:
                        result += "EXP\n";
                        break;
                }
                result += "Filling ticks:\t" + ((filled - start).Ticks / TimeSpan.TicksPerMillisecond).ToString("f1", CultureInfo.InvariantCulture) +
                    "\tSolving ticks:\t" + ((solved - filled).Ticks / TimeSpan.TicksPerMillisecond).ToString("f1", CultureInfo.InvariantCulture) +
                    "\tError:\t" + lengthError.ToString("f9", CultureInfo.InvariantCulture) + "\n";
                if (n <= 15) //if the number of points is less than 15, the solution and error vectors are added to the report
                {
                    result += "Original RHS Vector: " + VectorToString(rhs, n + m + 1, "f9") + "\n" +
                        "Solved RHS Vector: " + VectorToString(rhsSol, n + m + 1, "f9") + "\n" +
                        "Error RHS Vector: " + VectorToString(errorVector, n + m + 1, "f9") + "\n" +
                        "Solution of the LSE: " + VectorToString(sol, n + m + 1, "f9") + "\n";
                }
            }
            else         //Report for the analysis
            {
                result = "MODC:\t" + n.ToString() + "\t" + m.ToString() + "\t";
                switch (function)
                {
                    case RBFFunction.EXP:
                        result += "EXP\t";
                        break;
                    case RBFFunction.CP_C2:
                        result += "CP C2\t";
                        break;
                    case RBFFunction.TPS:
                        result += "TPS\t";
                        break;
                    case RBFFunction.MQ:
                        result += "MQ\t";
                        break;
                    default:
                        result += "EXP\t";
                        break;
                }
                result += ((filled - start).Ticks / TimeSpan.TicksPerMillisecond).ToString("f1", CultureInfo.InvariantCulture) + "\t" +
                    ((solved - filled).Ticks / TimeSpan.TicksPerMillisecond).ToString("f1", CultureInfo.InvariantCulture) +
                    "\t" + lengthError.ToString("f9", CultureInfo.InvariantCulture) + "\n";
            }
            #endregion
            return result;
        }

        /// <summary>
        /// Method for random creation of the test points.
        /// None of the points must not be equal to each other.
        /// </summary>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Array of created points.</returns>
        private double[][] CreatePoints(int n, int m)
        {
            double[][] result = new double[n][];
            int i, j;
            double xmin = -10.0, dx = 20.0;
            RandomN generator = new RandomN();

            for (i = 0; i < n; i++)
            {
                result[i] = new double[m];
                for (j = 0; j < m; j++)
                { result[i][j] = xmin + dx * generator.NextDouble();  }
            }
            return result;
        }

        /// <summary>
        /// Function to calculate interpolation point's value - multidimensional sphere.
        /// f(x)=sum_{i=1}^{n}(x_i)^2.
        /// </summary>
        /// <param name="x">Point.</param>
        /// <param name="m">Dimension of the point x.</param>
        /// <returns>Function value f(x).</returns>
        private double Function(double[] x, int m)
        {
            double result = 0.0;
            for (int i = 0; i < m; i++)
            { result += x[i] * x[i]; }
            return result;
        }
        #endregion

        #region Methods for solving the LSE
        /// <summary>
        /// LU decomposition and backsubstitution methods for solving LSE.
        /// </summary>
        /// <param name="matrix">Matrix of the LSE.</param>
        /// <param name="rhsVector">Right hand side vector of the LSE.</param>
        /// <param name="n">Number of equations in LSE.</param>
        /// <returns>Solution vector of the LSE.</returns>
        private double[] LUdecomposition(double[][] matrix, double[] rhsVector, int n)
        {
            double[] result = new double[n];

            #region Initialization
            double p = 0;
            double pom2;
            int k0 = 0, k, i, j;
            int pom1 = 0;
            double detOfP = 1;
            double[][] U = new double[n][];
            int[] pi = new int[n];
            double[][] L = new double[n][];
            double[] pom3;
            for (i = 0; i < n; i++)
            {
                L[i] = new double[n]; ;
                L[i][i] = 1.0;
                U[i] = new double[n];
                pi[i] = i;
            }
            U = matrix;
            #endregion

            #region Decomposition
            for (k = 0; k < n - 1; k++)
            {
                p = 0;
                for (i = k; i < n; i++)      // find the row with the biggest pivot
                {
                    if (Math.Abs(U[i][k]) > p)
                    {
                        p = Math.Abs(U[i][k]);
                        k0 = i;
                    }
                }

                pom1 = pi[k]; pi[k] = pi[k0]; pi[k0] = pom1;    // switch two rows in permutation matrix

                for (i = 0; i < k; i++)
                { pom2 = L[k][i]; L[k][i] = L[k0][i]; L[k0][i] = pom2; }

                if (k != k0) detOfP *= -1;

                pom3 = U[k]; U[k] = U[k0]; U[k0] = pom3;        // Switch rows in U

                for (i = k + 1; i < n; i++)
                {
                    L[i][k] = U[i][k] / U[k][k];
                    for (j = k; j < n; j++)
                    { U[i][j] = U[i][j] - L[i][k] * U[k][j]; }
                }
            }
            #endregion

            #region Permutation
            double[] b = new double[n];
            for (i = 0; i < n; i++)
            { b[i] = rhsVector[pi[i]];  }
            #endregion

            #region BackSubstitution
            double[] temp = new double[n];
            for (i = 0; i < n; i++)
            {
                temp[i] = b[i];
                for (j = 0; j < i; j++)
                { temp[i] -= L[i][j] * temp[j];  }
                temp[i] = temp[i] / L[i][i];
            }

            for (i = n - 1; i > -1; i--)
            {
                result[i] = temp[i];
                for (j = n - 1; j > i; j--)
                { result[i] -= U[i][j] * result[j];  }
                result[i] = result[i] / U[i][i];
            }
            #endregion

            return result;
        }

        /// <summary>
        /// Proposed decomposition method of the moment matrix of RPIM. Method also solves the LSE.
        /// It copies the solution in the result vector at the end. 
        /// </summary>
        /// <param name="Phi">Lower triangular part of the matrix Phi.</param>
        /// <param name="P">Matrix P.</param>
        /// <param name="f">Part of the right hand side vector of the LSE - vector f.</param>
        /// <param name="n">Dimension of the matrix Phi.</param>
        /// <param name="m">Number of columns of matrix P.</param>
        /// <returns></returns>
        private double[] modifiedCholeskyDecomposition(double[][] Phi, double[][] P, double[] f, int n, int m)
        {
            double[] result = new double[n + m + 1];

            #region Initialization
            //addition of the m due to the first monomial in the polynomial basis
            m++;

            int i, j, k, l;
            bool Aneg = false;
            double[] D = new double[n];
            double[][] A = new double[m][];
            double[] att = new double[n];
            double[] rhsL = new double[m];
            double[] b = new double[m];
            #endregion

            #region Cholesky decomposition of Phi and calculation of B. Matrix P becomes B in the process.
            D[0] = Phi[0][0];
            double sum = 0.0, sum2 = 0.0, temp = 0.0;
            double[] sumP = new double[m];
            for (i = 1; i < n; i++)
            { Phi[i][0] /= D[0]; }
            for (l = 0; l < m; l++)
            { P[l][0] /= D[0]; }
            for (i = 1; i < n; i++)
            {
                sum2 = Phi[i][0] * Phi[i][0] * D[0];
                for (l = 0; l < m; l++)
                { sumP[l] = P[l][0] * Phi[i][0] * D[0]; }
                for (j = 1; j < i; j++)
                {
                    sum = 0.0;
                    for (k = 0; k < j; k++)
                    { sum += D[k] * Phi[i][k] * Phi[j][k]; }
                    temp = Phi[i][j] - sum;
                    Phi[i][j] = (temp) / D[j];
                    sum2 += Phi[i][j] * temp;
                    for (l = 0; l < m; l++)
                    { sumP[l] += temp * P[l][j]; }
                }
                D[i] = Phi[i][i] - sum2;
                for (l = 0; l < m; l++)
                { P[l][i] = (P[l][i] - sumP[l]) / D[i]; }
            }
            #endregion

            #region Calculation of matrix A
            //*************************************//
            // Only the A[0][0] is calculated in   //
            // the frist step to determine the     //
            // definitnes of the A. If the A[0][0] //
            // is negative, then the A_{neg} is    //
            // calculated, which is later shown as //
            // Aneg=true. This effects in the      //
            // later calculation of att, b, etc.   //
            //*************************************//
            sum = P[0][0] * P[0][0];
            for (k = 1; k < n; k++)
            { sum += D[k] * P[0][k] * P[0][k]; }
            if (sum > 0.0)
            { Aneg = true; }
            if (Aneg)
            {
                for (i = 0; i < m; i++)
                {
                    A[i] = new double[m];
                    for (j = 0; j <= i; j++)
                    {
                        sum = P[i][0] * P[j][0] * D[0];
                        for (k = 1; k < n; k++)
                        { sum += D[k] * P[i][k] * P[j][k];  }
                        A[j][i] = sum;
                    }
                }
            }
            else
            {
                for (i = 0; i < m; i++)
                {
                    A[i] = new double[m];
                    for (j = 0; j <= i; j++)
                    {
                        sum = P[i][0] * P[j][0] * D[0];
                        for (k = 1; k < n; k++)
                        { sum += D[k] * P[i][k] * P[j][k]; }
                        A[j][i] -= sum;
                    }
                }
            }
            #endregion

            #region Calculation of the vector a'' - variable att[] and modified lower part of the rhs vector - vector 0.
            //Calculation of the att - forward substitution
            att[0] = f[0];
            for (i = 1; i < n; i++)
            {
                sum = f[i] - Phi[i][0] * att[0];
                for (k = 1; k < i; k++)
                { sum -= Phi[i][k] * att[k]; }
                att[i] = sum;
            }
            //Calculation of the modified lower part of the rhs vector - variable rhsL[]
            for (i = 0; i < m; i++)
            {
                sum = 0.0;
                for (j = 0; j < n; j++)
                { sum -= P[i][j] * att[j]; }
                rhsL[i] = sum;
            }
            #endregion

            #region Calculation of the part of solution vector b
            //Cholesky decomposition of matrix A
            for (i = 1; i < m; i++)
            { A[0][i] = A[0][i] / A[0][0]; }
            for (i = 1; i < m; i++)
            {
                sum2 = A[0][i] * A[0][i] * A[0][0];
                for (j = 1; j < i; j++)
                {
                    sum = 0.0;
                    for (k = 0; k < j; k++)
                    { sum += A[k][k] * A[k][i] * A[k][j]; }
                    temp = A[j][i] - sum;
                    A[j][i] = (temp) / A[j][j];
                    sum2 += A[j][i] * temp;
                }
                A[i][i] = A[i][i] - sum2;
            }
            //Forward substituion of the modified lower part of rhs vector - rhsL
            for (i = 1; i < m; i++)
            {
                sum = rhsL[0] * A[0][i];
                for (j = 1; j < i; j++)
                { sum += rhsL[j] * A[j][i]; }
                rhsL[i] -= sum;
            }
            //Backward substitution to calculate solution vector b
            for (i = m - 1; i > -1; i--)
            {
                b[i] = rhsL[i] / A[i][i];
                for (j = m - 1; j > i; j--)
                { b[i] -= b[j] * A[i][j]; }
            }
            #endregion

            #region Calculation of the vector at[] - it the att[]
            if (Aneg)
            {
                for (i = 0; i < n; i++)
                {
                    att[i] /= D[i];
                    sum = 0.0;
                    for (j = 0; j < m; j++)
                    { sum += P[j][i] * b[j]; }
                    att[i] += sum;
                }
            }
            else
            {
                for (i = 0; i < n; i++)
                {
                    att[i] /= D[i];
                    sum = 0.0;
                    for (j = 0; j < m; j++)
                    {  sum += P[j][i] * b[j]; }
                    att[i] -= sum;
                }
            }
            #endregion

            #region Calculation of the vector a[] (variable result[]) and copying vector b[] into the solution
            //Calculation of the vector a[]
            result[n - 1] = att[n - 1];

            for (i = n - 2; i > -1; i--)
            {
                sum = 0.0;
                for (j = n - 1; j > i; j--)
                { sum += Phi[j][i] * result[j]; }
                result[i] = att[i] - sum;
            }
            //Copying of the vector b[] into the solution vector
            if (Aneg)
            {
                for (j = 0, i = n; j < m; j++, i++)
                { result[i] = -b[j]; }
            }
            else
            {
                for (j = 0, i = n; j < m; j++, i++)
                { result[i] = b[j]; }
            }
            #endregion

            return result;
        }
        #endregion

        #region Methods for calculating whole matrix G0 - one method for each RBF function
        /// <summary>
        /// Method for calculating the G matrix in case of the EXP radial basis function.
        /// </summary>
        /// <param name="points">Array of the interpolationg points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated matrix G.</returns>
        private double[][] FillG_EXP(double[][] points, int n, int m)
        {
            double[][] G = new double[n + m + 1][];
            int i, j;
            for (i = 0; i < n; i++)
            {
                G[i] = new double[n + m + 1];
                for (j = 0; j <= i; j++)
                { G[i][j] = EXP(points[i], points[j], m); }
                for (j = i; j >= 0; j--)
                { G[j][i] = G[i][j]; }
            }
            for (i = n; i < n + m + 1; i++)
            { G[i] = new double[n + m + 1]; }
            for (i = 0; i < n; i++)
            {
                G[i][n] = 1.0;
                G[n][i] = 1.0;
                for (j = 0; j < m; j++)
                {
                    G[i][j + 1 + n] = points[i][j];
                    G[j + 1 + n][i] = points[i][j];
                }
            }
            return G;
        }

        /// <summary>
        /// Method for calculating the G matrix in case of the CP C2 radial basis function.
        /// </summary>
        /// <param name="points">Array of the interpolationg points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated matrix G.</returns>
        private double[][] FillG_CP_C2(double[][] points, int n, int m)
        {
            double[][] G = new double[n + m + 1][];
            int i, j;
            for (i = 0; i < n; i++)
            {
                G[i] = new double[n + m + 1];
                for (j = 0; j <= i; j++)
                { G[i][j] = CP_C2(points[i], points[j], m); }
                for (j = i; j >= 0; j--)
                { G[j][i] = G[i][j]; }
            }
            for (i = n; i < n + m + 1; i++)
            { G[i] = new double[n + m + 1]; }
            for (i = 0; i < n; i++)
            {
                G[i][n] = 1.0;
                G[n][i] = 1.0;
                for (j = 0; j < m; j++)
                {
                    G[i][j + 1 + n] = points[i][j];
                    G[j + 1 + n][i] = points[i][j];
                }
            }
            return G;
        }

        /// <summary>
        /// Method for calculating the G matrix in case of the TPS radial basis function.
        /// </summary>
        /// <param name="points">Array of the interpolationg points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated matrix G.</returns>
        private double[][] FillG_TPS(double[][] points, int n, int m)
        {
            double[][] G = new double[n + m + 1][];
            int i, j;
            for (i = 0; i < n; i++)
            {
                G[i] = new double[n + m + 1];
                for (j = 0; j <= i; j++)
                { G[i][j] = TPS(points[i], points[j], m); }
                for (j = i; j >= 0; j--)
                { G[j][i] = G[i][j]; }
            }
            for (i = n; i < n + m + 1; i++)
            { G[i] = new double[n + m + 1]; }
            for (i = 0; i < n; i++)
            {
                G[i][n] = 1.0;
                G[n][i] = 1.0;
                for (j = 0; j < m; j++)
                {
                    G[i][j + 1 + n] = points[i][j];
                    G[j + 1 + n][i] = points[i][j];
                }
            }
            return G;
        }

        /// <summary>
        /// Method for calculating the G matrix in case of the MQ radial basis function.
        /// </summary>
        /// <param name="points">Array of the interpolationg points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated matrix G.</returns>
        private double[][] FillG_MQ(double[][] points, int n, int m)
        {
            double[][] G = new double[n + m + 1][];
            int i, j;
            for (i = 0; i < n; i++)
            {
                G[i] = new double[n + m + 1];
                for (j = 0; j <= i; j++)
                { G[i][j] = MQ(points[i], points[j], m); }
                for (j = i; j >= 0; j--)
                { G[j][i] = G[i][j]; }
            }
            for (i = n; i < n + m + 1; i++)
            { G[i] = new double[n + m + 1]; }
            for (i = 0; i < n; i++)
            {
                G[i][n] = 1.0;
                G[n][i] = 1.0;
                for (j = 0; j < m; j++)
                {
                    G[i][j + 1 + n] = points[i][j];
                    G[j + 1 + n][i] = points[i][j];
                }
            }
            return G;
        }
        #endregion

        #region Methods for calculating lower triangular part of the Phi matrix - one for each RBF function
        /// <summary>
        /// Method for calculating the Phi matrix in case of the EXP radial basis function.
        /// </summary>
        /// <param name="points">Array of the interpolationg points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated lower triangular part of matrix PHI.</returns>
        private double[][] FillPHI_EXP(double[][] points, int n, int m)
        {
            double[][] G = new double[n][];
            int i, j;
            for (i = 0; i < n; i++)
            {
                G[i] = new double[n];
                for (j = 0; j <= i; j++)
                { G[i][j] = EXP(points[i], points[j], m);  }
            }
            return G;
        }

        /// <summary>
        /// Method for calculating the Phi matrix in case of the CP C2 radial basis function.
        /// </summary>
        /// <param name="points">Array of the interpolationg points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated lower triangular part of matrix PHI.</returns>
        private double[][] FillPHI_CP_C2(double[][] points, int n, int m)
        {
            double[][] G = new double[n][];
            int i, j;
            for (i = 0; i < n; i++)
            {
                G[i] = new double[n];
                for (j = 0; j <= i; j++)
                { G[i][j] = CP_C2(points[i], points[j], m); }
            }
            return G;
        }

        /// <summary>
        /// Method for calculating the Phi matrix in case of the TPS radial basis function.
        /// </summary>
        /// <param name="points">Array of the interpolationg points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated lower triangular part of matrix PHI.</returns>
        private double[][] FillPHI_TPS(double[][] points, int n, int m)
        {
            double[][] G = new double[n][];
            int i, j;
            for (i = 0; i < n; i++)
            {
                G[i] = new double[n];
                for (j = 0; j <= i; j++)
                { G[i][j] = TPS(points[i], points[j], m); }
            }
            return G;
        }

        /// <summary>
        /// Method for calculating the Phi matrix in case of the MQ radial basis function.
        /// </summary>
        /// <param name="points">Array of the interpolationg points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated lower triangular part of matrix PHI.</returns>
        private double[][] FillPHI_MQ(double[][] points, int n, int m)
        {
            double[][] G = new double[n][];
            int i, j;
            for (i = 0; i < n; i++)
            {
                G[i] = new double[n];
                for (j = 0; j <= i; j++)
                { G[i][j] = MQ(points[i], points[j], m); }
            }
            return G;
        }

        #endregion

        #region Methods for filling vectors and matricies that are same no matter which RBF function is used
        /// <summary>
        ///method for filling matrix P from the array of the interpolation points
        /// </summary>
        /// <param name="points">Array of m-dimensional points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Matrix P</returns>
        private double[][] FillPMatrix(double[][] points, int n, int m)
        {
            double[][] P = new double[m + 1][];
            int i, j, k;
            P[0] = new double[n];
            for (i = 0; i < n; i++)
            { P[0][i] = 1.0; }
            for (j = 0, k = 1; j < m; j++, k++)
            {
                P[k] = new double[n];
                for (i = 0; i < n; i++)
                { P[k][i] = points[i][j]; }
            }
            return P;
        }

        /// <summary>
        /// Function to construct right hand side vector of the LSE.
        /// Used in the defining the LSE for the LU decomposition and for the calculation of the errors.
        /// </summary>
        /// <param name="points">Array of points.</param>
        /// <param name="n">Number of points.</param>
        /// <param name="m">Dimension of points.</param>
        /// <returns>Calculated RHS vector.</returns>
        private double[] FillRHSVector(double[] points, int n, int m)
        {
            double[] RHS = new double[n + m + 1];
            int i = 0;
            for (i = 0; i < n; i++)
            { RHS[i] = points[i]; }
            for (i = n; i < n + m + 1; i++)
            { RHS[i] = 0.0; }
            return RHS;
        }

        /// <summary>
        /// Function to construct only f vector of the right hand side vector of the LSE.
        /// Used in the defining the LSE for the modified Cholesky decomposition.
        /// </summary>
        /// <param name="points">Array of points.</param>
        /// <param name="n">Number of points.</param>
        /// <returns>Calculated f vector.</returns>
        private double[] Fill_fVector(double[] points, int n)
        {
            double[] f = new double[n];
            int i = 0;
            for (i = 0; i < n; i++)
            { f[i] = points[i]; }
            return f;
        }
        #endregion

        #region RBF Functions
        //Parameters of the RBF functions
        private double CP_C2_rr = 9.0;
        private double MQ_c = 0.8;
        private double MQ_aa = 0.3;
        private double EXP_c = 0.1;

        /// <summary>
        /// TPS radial basis function.
        /// </summary>
        /// <param name="point1">Point i.</param>
        /// <param name="point2">Point j.</param>
        /// <param name="m">Dimension of the points.</param>
        /// <returns>Value of the TPS RBF function.</returns>
        private double TPS(double[] point1, double[] point2, int m)
        {
            double result = 0.0, r = 0.0;
            for (int i = 0; i < m; i++)
            { r += (point1[i] - point2[i]) * (point1[i] - point2[i]); }
            if (r != 0.0)
            { result = r * Math.Log(Math.Sqrt(r))+1.0; }
            else { result = 1.0; }
            return result;
        }

        /// <summary>
        /// EXP radial basis function.
        /// </summary>
        /// <param name="point1">Point i.</param>
        /// <param name="point2">Point j.</param>
        /// <param name="m">Dimension of the points.</param>
        /// <returns>Value of the EXP RBF function.</returns>
        private double EXP(double[] point1, double[] point2, int m)
        {
            double result = 0.0, r = 0.0;
            for (int i = 0; i < m; i++)
            { r += (point1[i] - point2[i]) * (point1[i] - point2[i]); }
            result = Math.Exp(-EXP_c * r);
            return result;
        }

        /// <summary>
        /// MQ radial basis function.
        /// </summary>
        /// <param name="point1">Point i.</param>
        /// <param name="point2">Point j.</param>
        /// <param name="m">Dimension of the points.</param>
        /// <returns>Value of the MQ RBF function.</returns>
        private double MQ(double[] point1, double[] point2, int m)
        {
            double result = 0.0, r = 0.0;
            for (int i = 0; i < m; i++)
            { r += (point1[i] - point2[i]) * (point1[i] - point2[i]); }
            result = Math.Pow((r + MQ_aa), MQ_c);
            return result;
        }

        /// <summary>
        /// CP_C2 radial basis function.
        /// </summary>
        /// <param name="point1">Point i.</param>
        /// <param name="point2">Point j.</param>
        /// <param name="m">Dimension of the points.</param>
        /// <returns>Value of the CP_C2 RBF function.</returns>
        private double CP_C2(double[] point1, double[] point2, int m)
        {
            double result = 0.0, r = 0.0;
            for (int i = 0; i < m; i++)
            { r += (point1[i] - point2[i]) * (point1[i] - point2[i]); }
            r = Math.Sqrt(r);
            if (r < CP_C2_rr)
            {
                r /= CP_C2_rr;
                result = Math.Pow((1 - r), 4.0) * (4.0 * r + 1);
            }
            else { result = 0.0; }
            return result;
        }
        #endregion

        #region Methods for error calculation
        /// <summary>
        /// Method for calculation of the error vector
        /// </summary>
        /// <param name="m">Matrix of the original LSE.</param>
        /// <param name="sol">Numerical solution of the LSE.</param>
        /// <param name="n">Size of the LSE.</param>
        /// <param name="rhs">Exact right hand side vector of the LSE.</param>
        /// <returns>The vector containig the errors for each component.</returns>
        private double[] CalculateRHSVector(double[][] m, double[] sol, int n)
        {
            double[] result = new double[n];
            int i, j;
            for (i = 0; i < n; i++)
            {
                result[i] = 0.0;
                for (j = 0; j < n; j++)
                { result[i] += m[i][j] * sol[j]; }
            }
            return result;
        }

        /// <summary>
        /// Function to calculate the difference between the original and the solved RHS vectors.
        /// </summary>
        /// <param name="solution">RHS vector obtained by multiplying original matrix and the solution vector.</param>
        /// <param name="original">RHS vector used in the calculation of the solution.</param>
        /// <param name="n">Lengths fo the vectors.</param>
        /// <returns>Vector containing the differences of vectors.</returns>
        private double[] CalculateErrorVector(double[] solution, double[] original, int n)
        {
            double[] result = new double[n];
            for (int i = 0; i < n; i++)
            { result[i] = original[i] - solution[i]; }
            return result;
        }

        /// <summary>
        /// Method for calculation of the lenght of the vector.
        /// </summary>
        /// <param name="vec">Vector.</param>
        /// <param name="n">Dimension of the vector.</param>
        /// <returns>The lenght of the vector vec.</returns>
        private double LengthVector(double[] vec, int n)
        {
            double result = 0.0;
            for (int i = 0; i < n; i++)
            { result += vec[i] * vec[i]; }
            result = Math.Sqrt(result);
            return result;
        }

        /// <summary>
        /// Returns the string representation of the vector.
        /// </summary>
        /// <param name="vec">Vector to convert to string.</param>
        /// <param name="n">Lenght of the vector.</param>
        /// <param name="format">Number format for the conversion.</param>
        /// <returns>String representation of the vector.</returns>
        private string VectorToString(double[] vec, int n, string format)
        {
            string result = "{";
            int i;
            for (i = 0; i < n - 1; i++)
            { result += vec[i].ToString(format, CultureInfo.InvariantCulture) + ", "; }
            result += vec[i].ToString(format, CultureInfo.InvariantCulture) + "}";
            return result;
        }
        #endregion
    }

    #region Object for creatign random numbers between 0 and 1
    internal class RandomN
    {
        public RandomN() { }

        public double NextDouble()
        {
            byte[] random_bytes = new byte[8];
            RNGCryptoServiceProvider rng = new RNGCryptoServiceProvider();
            rng.GetNonZeroBytes(random_bytes);
            double result = 0.0;
            try
            { result = Math.Abs(((double)BitConverter.ToUInt32(random_bytes, 0)) / (1.0 + (double)UInt32.MaxValue)); }
            catch
            {
                new RNGCryptoServiceProvider().GetBytes(random_bytes);
                result = Math.Abs(((double)BitConverter.ToUInt32(random_bytes, 0)) / (1.0 + (double)UInt32.MaxValue));
            }
            return result;
        }
    }
    #endregion
}
